https://www.neonscience.org/resources/learning-hub/tutorials/phenocam-phenor-modeling
For the latest version, install phenor from here: https://github.com/bluegreen-labs/phenor (be sure to
update all packages for installation to work, and install GenSA)
# example from website: Bartlett NH
# phenocamr::download_phenocam(
# frequency = 3,
# veg_type = "DB", #DB for decid broadlead, EN for evergreen, default is ALL
# roi_id = 1000,
# site = "bartlettir",
# phenophase = TRUE,
# out_dir = "." # saves in current working dir
# )
# evergreens in del norte county (close to north humboldt)
# phenocamr::download_phenocam(
# frequency = 3,
# veg_type = "EN", #DB for decid broadlead, EN for evergreen, default is ALL
# roi_id = 1000,
# site = "delnortecounty1",
# phenophase = TRUE,
# out_dir = "." # saves in current working dir
# )
# shrubs in crownpoint NM, NE of Gallup at Navajo Technical University
#phenocamr::download_phenocam(
# frequency = 3,
# roi_id = 1000,
# site = "ntupinyongarden",
# phenophase = TRUE,
# out_dir = "." # saves in current working dir
# )
# two files: time series and transition data
Find the metadata here: https://phenocam.nau.edu/webcam/tools/summary_file_format/
https://bluegreen-labs.github.io/phenocamr/articles/phenocamr-vignette.html
source has clear steps, but uses base plotting
# read as phenocamr type
df = read_phenocam("example_data/bartlettir_DB_1000_3day.csv")
df <- expand_phenocam(df) # expand every 3 day to every day
df <- detect_outliers(df) # detect outliers
df <- smooth_ts(df) # smooth data using AIC
phenology_dates <- phenophases(df, internal = TRUE) # calculate rising and falling
as.data.frame(df$data) %>%
mutate(date = as.Date(date)) %>%
ggplot(aes(x=date, y=smooth_gcc_90)) +
geom_line() +
labs(x="Date", y="GCC") +
# rising "spring" greenup dates
geom_vline(data=phenology_dates$rising, aes(xintercept=transition_50), col="green") +
# falling 'autumn' senescence dates
geom_vline(data=phenology_dates$falling, aes(xintercept=transition_50), col="brown")
Use the function merge_daymet
Daymet is continuous, gridded estimates of weather and climate variables creted through statistical interpolation of ground based measurements. Read more here: https://daymet.ornl.gov/
bnh_clim = merge_daymet("example_data/bartlettir_DB_1000_3day.csv")
names(bnh_clim)
## [1] "site" "veg_type" "roi_id" "frequency"
## [5] "lat" "lon" "elev" "solar_elev_min"
## [9] "header" "data"
bnh_data = bnh_clim$data
# plot daily rainfall and temperatures
ggplot(bnh_data, aes(x=date, y=prcp..mm.day.))+
geom_line() +
labs(y="precip [mm]", title="Daily Precipitation") +
theme_bw()
ggplot(bnh_data) +
geom_point(aes(x=date, y=tmax..deg.c., col="tmax"), col="red") +
geom_point(aes(x=date, y=tmin..deg.c., col="tmin"), col="blue") +
labs(y="temp. deg C", title = "Daily temperature") +
theme_bw()
# day of the year patterns
ggplot(bnh_data) +
geom_point(aes(x=doy, y=gcc_90))
## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_point()`).
# precipitation
ggplot(bnh_data) +
geom_point(aes(x=prcp..mm.day., y=gcc_90, col=doy)) + scale_color_viridis_c()
## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_point()`).
# day of the year with daylength
ggplot(bnh_data) +
geom_point(aes(x=doy, y=gcc_90, col=dayl..s.)) + scale_color_viridis_c()
## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_point()`).
# daylength
ggplot(bnh_data) +
geom_point(aes(x=doy, y=gcc_90, col=dayl..s.)) + scale_color_viridis_c()
## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_point()`).
# min temperatures
ggplot(bnh_data) +
geom_point(aes(x=tmin..deg.c., y=gcc_90))
## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_point()`).
# max temperatures
ggplot(bnh_data) +
geom_point(aes(x=tmax..deg.c., y=gcc_90))
## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_point()`).
# subset data
data_mod = bnh_data %>% dplyr::filter(year<2014)
# start with one variable
lm_test = lm(data=data_mod, gcc_90~dayl..s.)
# see residuals, p-vals, Rsquares, etc.
summary(lm_test)
##
## Call:
## lm(formula = gcc_90 ~ dayl..s., data = data_mod)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.072102 -0.011738 0.008489 0.015394 0.053883
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.042e-01 4.963e-03 41.15 <2e-16 ***
## dayl..s. 4.627e-06 1.123e-07 41.21 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02429 on 672 degrees of freedom
## (22 observations deleted due to missingness)
## Multiple R-squared: 0.7165, Adjusted R-squared: 0.7161
## F-statistic: 1699 on 1 and 672 DF, p-value: < 2.2e-16
# check with more variables
lm_test2 = lm(data=data_mod, gcc_90~dayl..s.+tmin..deg.c.)
# predict with remaining years
data_pred = bnh_data %>%
dplyr::filter(year>=2014) %>%
select(gcc_90, dayl..s., tmin..deg.c., year, doy, date)
# predict with one var
data_pred$pred1_gcc_90 = predict.lm(lm_test, data_pred)
# predict with 2 vars
data_pred$pred2_gcc_90 = predict.lm(lm_test2, data_pred)
ggplot(data_pred) +
geom_point(aes(x=date, y=pred1_gcc_90, col="1 var")) +
geom_point(aes(x=date, y=pred2_gcc_90, col="2 vars"))+
geom_point(aes(x=date, y=gcc_90, col="obs"))
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
Transition dates are 25% of GCC amplitude - included in phenocamr data structure
Ti = mean daily temps
Li = daylength
# read in parameter ranges
path <- sprintf("%s/extdata/parameter_ranges.csv",path.package("phenor"))
par_ranges <- read.table(path,
header = TRUE,
sep = ",")
lin_par = par_ranges %>%
dplyr::filter(model=="LIN") %>%
select(a,b)
# load data from phenor package
bnh = phenocam_DB$bartlettir # in proper format for models
# plot temps
tmin = data.frame(bnh$Ti)
colnames(tmin) <- bnh$year
tmin$doy = bnh$doy
# shifting doy so that 0 is beginning of year jan 1
tmin$doy = ifelse(tmin$doy<0,
tmin$doy+365,
tmin$doy)
tmin2 = pivot_longer(tmin, names_to="year", values_to="tmin", cols=-doy)
# plot temperatures
tmin2 %>%
ggplot(aes(x=doy, y=tmin, col=year)) +
geom_line()
# make year / spring dates
spring = data.frame(year=bnh$year, greenup=bnh$transition_dates)
# plot this pattern
ggplot(spring, aes(x=year, y=greenup)) +
geom_point() +
geom_smooth(method="lm")
## `geom_smooth()` using formula = 'y ~ x'
try_lm = lm(data=spring, greenup~year)
# add dates to plot
tmin2 %>%
ggplot(aes(x=doy, y=tmin, col=as.factor(year))) +
geom_line() +
geom_vline(data=spring, aes(xintercept=greenup, col=as.factor(year))) + scale_color_viridis_d()
# plot spring df with tmin
ann_tmin = tmin2 %>% group_by(year) %>%
summarise(tmin_min = min(tmin),
tmin_max = max(tmin),
tmin_mn = mean(tmin)) %>%
mutate(year=as.double(year)) %>%
inner_join(spring, by="year")
ggplot(ann_tmin, aes(x=year,
y=tmin_mn,
ymin=tmin_min,
ymax=tmin_max)) + geom_pointrange()
ggplot(ann_tmin, aes(y=greenup,
x=tmin_mn)) +
geom_point() + geom_smooth(method='lm')
## `geom_smooth()` using formula = 'y ~ x'
# could also look at spring temps
# we will define spring as 60-151 (March-Apr-May)
ann_tmin_spr = tmin2 %>%
filter(doy>=60 &
doy<=151) %>%
group_by(year) %>%
summarise(tmin_min = min(tmin),
tmin_max = max(tmin),
tmin_mn = mean(tmin)) %>%
mutate(year=as.double(year)) %>%
inner_join(spring, by="year")
# plot with mean min tmps
ggplot(ann_tmin_spr, aes(y=greenup,
x=tmin_mn)) +
geom_point() + geom_smooth(method='lm')
## `geom_smooth()` using formula = 'y ~ x'
# estimate params by working backwards with
spr_lm = lm(data=ann_tmin_spr, greenup~tmin_mn)
summary(spr_lm)
##
## Call:
## lm(formula = greenup ~ tmin_mn, data = ann_tmin_spr)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.158 -1.864 -1.031 2.016 3.275
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 138.8544 3.6532 38.009 2.21e-08 ***
## tmin_mn -2.3237 0.6165 -3.769 0.0093 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.517 on 6 degrees of freedom
## Multiple R-squared: 0.7031, Adjusted R-squared: 0.6536
## F-statistic: 14.21 on 1 and 6 DF, p-value: 0.009299
# data must be formatted for input
# LIN runs a linear model of spring mean temps
# doy = a*mean_spring_temp + b
tmp = LIN(par=c(a=-2.3, b=138.8), data=bnh)
ann_tmin_spr$lm_pred = tmp
ggplot(ann_tmin_spr) +
geom_point(aes(x=greenup, y=lm_pred)) +
geom_abline(aes(intercept=0, slope=1)) +
labs(x="observed", y="predicted")